Non-BCS superfluidity in trapped ultracold Fermi gases 
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Superconductivity and superfluidity of fermions require, within the BCS theory, matching of the 
Fermi energies of the two interacting Fermion species. Difference in the number densities of the 
two species leads either to a normal state, to phase separation, or - potentially - to exotic forms of 
superfluidity such as FFLO-state, Sarma state or breached pair state. We consider ultracold Fermi 
gases with polarization, i.e. spin-density imbalance. We show that, due to the gases being trapped 
and isolated from the environment in terms of particle exchange, exotic forms of superfluidity appear 
as a shell around the BCS-superfluid core of the gas and, for large density imbalance, in the core 
as well. We obtain these results by describing the effect of the trapping potential by using the 
Bogoliubov-de Gennes equations. For comparison to experiments, we calculate also the condensate 
fraction, and show that, in the center of the trap, a polarized superfluid leads to a small dip in the 
central density difference. We compare the results to those given by local density approximation 
and flnd qualitatively different behavior. 

PACS numbers: 03.75.Ss, 03.75.Hh, 05.30.Jp 



INTRODUCTION 



There are several suggestions of non-BCS superfluidity 
for fermion systems with spin-population imbalance, i.e. 
with the polarization P = {N^-Ni)/{Ni+Ni) ^ where 
Na are the particle numbers of the (pseudo) spins. The 
FFLO-state [l|, [ah the Sarma state Q and the breached 
pair (BP) state [J| all appear as extremal points of the 
mean-field energy of the system. Such superfiuidity is 
of interest e.g. in condensed matter, high-energy and 
nuclear physics [5!], but firm experimental evidence is 
lacking. With the recently realized superfiuids of alkali 
Fermi gases, the first studies of density- imbalanced gases 
[6, 7, 8, 9] have shown that these systems offer unprece- 
dented opportunities for investigating this question. The 
FFLO-state is predicted to appear in a narrow parame- 
ter window in several systems |lOl . I 111 . Il2l . Il3l |. and the 
Sarma/BP state has been shown to be unstable under 
many conditions |l4l . Il5| . The existence of these exotic 
forms of superfiuidity is thus an intriguingly subtle ques- 
tion and requires careful analysis, taking into account 
the specific features of the physical system. In this arti- 
cle, we demonstrate the important consequences of such 
features in case of trapped ultracold alkali gases. First, 
in ultracold gases, the physical system under study is 
isolated from the environment in terms of particle ex- 
change. This could be contrasted to an electronic system 
where external voltage fixes the chemical potential by 
allowing particle exchange between the system of inter- 
est and the environment. Second, the gas is trapped by 
an inhomogeneous potential (often of harmonic form). 
The second issue leads to phase separation of super- 
fiuid and normal phases, as shown by a series of ex- 




and theoretical studies 
The finiteness and 



periments 

harmonic confinement of the system leads, as we show in 
this article, to stabilization of exotic forms of superfiu- 
idity in a shell surrounding a BCS-superfiuid core, and 



eventually in the core itself. 

Ultracold Fermi gases offer the possibility of realizing 
the BCS-BEC crossover with the use of Feshbach reso- 
nances. At the resonance point, the system is a strongly 
interacting Fermi superfiuid, and on different sides away 
from the resonance it is either a BEC of molecules or 
a weakly interacting Fermi gas realizing a BCS super- 
fiuid. Condensation of molecules, fermion pairs, pairing 
gap and the crossover behavior have been studied by a se- 
ries of experiments and vortices confirming superfiuidity 
have been created, for a review see for instance [23|. Re- 
cently, seminal experiments studying density imbalanced 
(nonzero polarization P) Fermi gase s were done [3, 0, S]- 
The analysis presented in [a, [la, |l7| combined the study 
of vortex patterns, density profiles (in [l7| 3D recon- 
structed) and condensate fractions, and thereby showed 
that, in the trapped system, a superfiuid core (with equal 
densities of the components t and J,) appears in the cen- 
ter of the trap and the excess atoms of the majority 
component (t) tend to be located on the edges of the 
trap. In the following, the BCS paired superfiuid (SF) 
implies equal local densities, whereas the polarized su- 
perfiuid (PS) at zero temperature denotes exotic forms 
of superfiuidity with non-zero local density difference of 
the t and | components. In addition, the normal state 
(N) shell surrounding the condensate can either be par- 
tially polarized or fully polarized depending on whether 
the minor component is present or not. 

In connection to the experiments p, 0, S| , several au- 
thors have analyzed the trapped, polarized Fermi gas us- 
ing the local density approximation (within mean-field 
theory) [13, M, M, M, Hi, M, M,m^- In these 
works, the problem was solved at each spatial point, ap- 
plying a local chemical potential, and the stability of the 
solutions at that point was determined by their grand 
potential energies. This leads, at zero temperature, to 
the exclusion of Sarma/BP state which is known to be a 
maximum point of the grand potential for fixed chemical 
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Figure 1: (color online) The gap A(r) and density difference 
Sn = n^{r) - n|(r) profiles for the LDA and BdG calculations. 



Figure 2: (color online) The density profiles of the n|(r) for 
BdG (solid) and LDA (dashed), and ni{r) for BdG (dot- 
dashed) and LDA (dotted). 



potentials. This gives the following qualitative picture: a 
superfluid (SF) core, with equal densities, appears at the 
center of the trap, and is surrounded by a normal state 
(N) shell (only on the BEC side of the Feshbach resonance 
a coexistence of molecular BEC and free atoms was seen, 
as obviously expected since molecule formation does not 
require matching Fermi energies). This is, however, a 
rather approximative description of the system because 
LDA assumes a smooth variation of the density. Since the 
studies [l3,|23,[23,[2i,[2l,[2i| imply a step in the density 
at the interface of the SF core and the N shell one can 
question the validity of LDA at that boundary, as also re- 



undary, as also re- 

210, Q, Sin 



marked by many of the authors of] 

and recent beyond-LDA studies |l8l . I32l . I3S 

The existence of a PS for trapped, strongly inter- 
acting gases was indicated already in our earlier study 
|2l| . where the treatment of the trapping potential by 
BogoHubov-de Gennes (BdG) equations revealed oscil- 
lations of the order parameter in an area located be- 
tween the SF core and the N shell. Such oscillations re- 
semble the nonuniform order parameter associated with 
the FFLO state, therefore we refer to these oscillations 
as "FFLO-type state". In this article, we show that 
the BdG analysis predicts such polarized superfluid in 
trapped Fermi gases not only as a shell effect but, for 
large polarizations, as a feature that extends through the 
whole system. We confirm that superfiuidity and finite 
local density difference indeed co-exist in the center of 
the trap by calculating the condensate fraction, central 
gap and the core polarization. In previous works using 
BdG for trapped gases [ill, [l3, UM, H] such an anal- 
ysis, confirming the existence of a polarized superfluid 
in the center, was not performed. Moreover, the works 
[ill . [l2 . [34 | considered the weak interaction limit whereas 
we have extended the BdG calculation to the unitarity 
regime and actually consider the whole crossover from 



the BCS to BEC side. We analyze the dependence of 
the oscillations on the system size (particle number) and 
discuss the connection to phenomena occurring in super- 
conductor - ferromagnet interfaces. For comparison, we 
also perform calculations using LDA. The LDA analysis 
predicts a polarized superfluid only at finite temperature, 
not at zero temperature like the BdG calculation. There- 
fore our results show that, in trapped Fermi gases, LDA 
has to be applied with care; beyond LDA approaches are 
not only needed for describing the interfaces of the sys- 
tem properly but, in some cases, features not predicted 
by LDA may become dominant. 

This article is structured as follows. In section II. we 
present a review of the BogoHubov - de Gennes (BdG) ap- 
proach for describing pairing at the mean field level. We 
discuss the use of local density approximation, and the 
expansion in harmonic oscillator states as special cases of 
this general scheme. In the rest of the paper, we use the 
following terminology: BdG refers to the case where har- 
monic oscillator state basis has been used, and LDA to 
the use of local density approximation. The results are 
presented in section II. The appearance of the FFLO- 
type oscillations is discussed in subsection II. A and the 
case of large polarizations, when such oscillations span 
the whole system, is discussed in section II. B. We also 
analyze the dependence of these features on the interac- 
tion strength through the BCS-BEC crossover (section 
II. A. 1.) and on the system size (atom number) (section 
II. A. 2). The condensate fraction is calculated in section 
II. C. and the contributions of different harmonic oscil- 
lator states to the pairing are discussed in section II. D. 
Conclusions and discussion are presented in section III. 
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Figure 3: (color online) The gap and polarization profiles at 
slightly elevated temperature T/Tp — 0.05 for the LDA calcu- 
lations and at low polarization. In contrast to the zero tem- 
perature case the thermal excitation in the BCS core starts 
to show up as a tail of the density difference into the core. 



II. BOGOLIUBOV - DE GENNES APPROACH 

In order to properly account for the inhomogeneity due 
to the presence of the trap potential we will in the follow- 
ing use the Bogoliubov-de Gennes equations [35|, which 
in a more general context are also known as the self- 
consistent Hartree-Fock-Bogoliubov equations. The im- 
balanced two-component Fermi gas is described by the 
grand canonical Hamiltonian 



H = Y,jd^r¥,{r) 



2y72 



2m 



/V + y(r) *,(r) 



d^rifr' *j(r)^|(r)t/(r-r')^i(r)^t(r), (1) 



where *iE'o-(r), ^J.(r) are the real space annihilation and 
creation operators for an atom with spin a at position 
r, ficr is the chemical potential for the components a, 
V{r) = ^muj'^r'^ is the external (isotropic and harmonic) 
trapping potential, and U{r) — US{r) is the interatomic 
atom-atom contact interaction potential. We apply the 
contact potential interaction and the mean-field approx- 
imation 



C/*|(r)*j(r)l'x(r)#t(r) 



|A(r) 
U 



C/?i|(r)n|(r) 



+ C/n|(r)4'J(r)*|(r) + L/nx(r)*|(r)4'i-(r) 

+ A(r)*|(r)*|(r) + A*(r)*|(r)*|(r), (2) 

where A(r) = (^|(r)*|(r)) and n^{r) = (#t.(r)#^(r)). 
The first two terms are simply numbers, so they are ne- 
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J2 J d^'r Un^{r)i'Ur)i'a{r) 
(fr(A{r)i/\{r)ii\{r)+ H. 



(3) 



with a being the other component of a. The second line 
is the mean-field Hartree corrections which represents 
an asymmetric shift to the chemical potential ^a of the 
gas of a atoms which is proportional to the density of 
the other component Ug-. In the Bogoliubov-de Gennes 
approach one expands the field operators on an appro- 
priate basis set, dictated by the symmetries of the non- 
interacting part of the Hamiltonian and usually charac- 
terized by a set {ij} of good quantum numbers, in order 
to diagonalize the Hamiltonian as in the uniform case. 
In the balanced non-uniform case and in the absence of 
superfiow the generaHzed BCS pairing is in general be- 
tween atoms in time-reversed states. The problem can 
now be solved by introducing the generalized canonical 
BogoHubov-Valatin transformation to the new fermionic 
operators a^,ajj which amounts to the expansion 



where the overhead bar designates the quantum numbers 
of the time-reversed state for r] and the other hyperfine 
spin state for cr, and with s-f = l,si = —1. We note 
that in the analogous case of electronic pairing in super- 
conductors a denote the time-reversed spin part of the 
wavefunction. Prom the requirement that the new oper- 
ators arj,alj diagonalize the Hamiltonian ^ one derives 
the matrix Bogoliubov-de Gennes equation [35|| 



{Hofs - Afi)y),, = £',,<p,,(r). 



(4) 



for the spinor (p{r) = (u,,(r), u^(r)), where 77 denotes 
a set of appropriate quantum numbers and Urj,Vjj are 
therefore to be regarded as subspinors, Hq is the non- 
interacting diagonal part of the Hamiltonian, potentially 
including a trapping potential and Hartree shifts to the 
chemical potentials, A is the pairing field part of the 
Hamiltonian, Er, is the eigenenergy. The products with 
the Pauli matrices Ti on the left hand side of Eq. ^ are 
to be understood as a direct products. The selfconsistent 
aspect of the method is due to the fact that the chemical 
potentials, mean-field Hartree and pairing fields are to be 
selfconsistently determined through the gap and number 
equations. The details of the BdG calculation for the 
harmonic trap eigenstates are presented in Section IIIBI 
Next we turn to the discussion of the local density ap- 
proximation. 



A. Local density approximation 

We first consider the local density approximation 
which is assumed to be valid for sufficiently large con- 
densates. The starting point for the LDA calculation is 
to solve the problem in the uniform case. The transla- 
tional invariance of a uniform superffuid impHes that the 
plane wave states ^(r) = V~^/^ ^^ e^^"^ak can be used 
to diagonalize the Hamiltonian. The main assumption of 
the LDA is that the system is locally homogeneous and 
therefore we initially consider the Hamiltonian H m ^ 
for a uniform system (i.e. V{r) = 0) with contact inter- 
actions, which in the momentum representation reads 



H = >,eka\^ako 



U 



\^J J 



k.a 



k,k' 



(5) 
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where aj,^ , Uka are the creation and annihilation oper- 
ators for free atoms with momentum hk, and the ki- 
netic energy e^ = h^k^/{2m). The chemical potentials 
fia- are introduced as Lagrange multipliers for the par- 
ticle numbers A^^r- We define the average chemical po- 
tential /i — (ni + /i|)/2 and the imbalance potential 
Sfj, = {fi^ — fii)/2, such that fi^ = fJ. + s^S^j,, with Sa 
defined as above. Within the mean field approxima- 
tion the Hamiltonian can be diagonaHzed by the stan- 
dard Bogoliubov-Valatin transformation. In the uniform 
case the order parameter is the usual A = {uk^a^ki) 
which satisfies the gap equation. Within the semiclas- 
sical approximation the gap A(r) is assumed to depend 
weakly on the radial position r — \r\ and, as is com- 
monly done, we introduce the local chemical potential 
Mo- ~^ Mcr('') — ^'a — y{''')i where V{r) is the trapping 
potential. The quasi-particle dispersions depend on the 
gap and therefore on the radial position through the local 
chemical potentials and the gap profile. In the most gen- 
eral case the quasi-particle dispersion relation becomes 

Eka{r) - \ {ika{r) - ika{r)) + Ekir), 

where Ek{r) = [^k{r) + A^r)]^/^, with ^kir) ^ 
i^kAr) + ^ksir)) /2, and^kair) = h^k^/{2m)-ii,{r). As 
we will assume pairing between atoms with equal mass we 
get Eka = -Scr6fi + Ek, and therefore 2Ek = Eka + Eks, 
and ^fc = Sk — fJ-ir) with ^, 5^ defined above. Within 
LDA, both the gap and the local chemical potentials 
depend on the radial position and therefore the system 
may in general be locally polarized. We define an aver- 
aged Fermi energy scale E-p from the total atom num- 
ber N = [Ey / {huj)]^ / i , for a balanced non-interacting 
Fermi gas in a harmonic trap potential V(r) = TOW^r^/2, 
where w is the trap frequency. The characteristic scale 
of the size of the cloud is the Thomas-Fermi radius is 
i?TF = [2£;F/("itj2)]i/2_ Within LDA the local gap equa- 
tion at position / reads 



1 = 



V 



E 



1 I - nFJEk^jr)) - nyjEkiir)) 

2£k 2Ek{r) 



with C/* = Airh'^as/m being the usual regularized effec- 
tive interaction strength which replaces the bare interac- 
tion U, and Qs is the s-wave scattering length, and where 
we have regularized the ultraviolet divergence appearing 
momentum integrals arising from unphysical properties 
of the contact interaction [3a|- The gap equation is an 
implicit equation for the gap profile A(r) [3^ which for 
an isotropic trap is only a function of the radial position 
r — \r\. The number equation and the density profiles 
are determined from the thermodynamic relation 



iV„ 






(fir n^{r), 



where the radial distribution for the a component is 

"''^('') " vE ["^liT-FiEkair)) +vlTiF{-Eks{r))] , 



V 



with a denotes the other component of a. The polariza- 
tion P is 



P = 






1 

N 



d r 6n{r), 



(7) 



(6) 



where Sn{r) = n^{r) — nj^(r). The equations for the po- 
larization P and for the minority component N^/N = 
(1 — P)/2 are numerically solved by iteration together 
with the gap profile as follows. We first make an ini- 
tial guess for the values of the chemical potentials /i and 
Sfi at fixed coupling kpas- The r dependent gap A(r), 
minor density ni{r) and polarization density Sn{r) are 
then discretized on a radial grid of sufficient resolution. 
From the initial values of n and 5^ and at fixed cou- 
pling the gap profile A(r) is calculated by solving the 
gap equation |[6]) as a root problem. The gap profile is 
then subsequently used to calculate the densities ni{r) 
and Sn{r) as functions of /i and (5/i. The minor number 
and and polarization equations are then solved as a mul- 
tidimensional root problem for the average chemical po- 
tential /i and the bare depairing width S^. The system of 
equations are iterated until the gap profile and the chem- 
ical potentials are sufficiently converged and then finally 
n^(r) — Sn{r) +ni{r). The present method applied here 
is well suited for considering effects of the Hartree terms, 
but we have not included those in the results presented 
here. 



B. Harmonic trap eigenstates 

The Bogoliubov-de Gennes (BdG) approach allows 
treating the effect of the harmonic trapping potential ex- 
actly. The approach has been used for polarized Fermi 
gases [11, 12, 21, 34] because it avoids some of the prob- 
lems of a local density approximation approach. We have 
used it in our previous work [2l| and give in this section 
a detailed description of the calculations on which the 
results in [2l| and in this article are based. One particu- 
larly important advantage of this approach is the proper 



description of interface effects in a phase separated gas, 
following from the nonlocal nature of the solution and the 
wavefunctions. Eventually, in the limit of a large num- 
ber of atoms, the LDA and BdG solutions are expected 
to become the same. 

The BdG approach used here is a generaHzation, to 
the imbalanced densities case, of the calculation outHned 
in Ref. [38] which was for an unpolarized Fermi gas. 
For the imbalanced case the generalization amounts to 
introducing different chemical potentials ^x and ^| for 
the two species of fermionic atoms. We now expand the 
wavefunctions ^(j(r) in the basis of the 3D harmonic 
oscillator eigenstates 



and the pairing field is described by 



^a{r) 



E 



Rnlir)Yim{r)anlma, 



(8) 



where the quantum numbers {r]} = {n, I, m} are the ra- 
dial quantum number n counting the nodes in the radial 
function, I is the orbital angular momentum and m is 
the projected angular momentum onto the axis of quan- 
tization. Here Yim{r) are the spherical harmonics with 
f = {6, if) and the radial part of the wavefunction is 



Rni{r) 



(9) 

where r — r/uosc with the harmonic oscillator length 

being Oosc = [f^/ifn^)]^^^, and L„ (f^) is an asso- 
ciated Laguerre polynomial. The characteristic scale 
of the cloud is given by the Thomas-Fermi radius 
i?|p = 2EF/{moj^) and therefore Rtf = {2AN)^^^aosc. 
The spherical symmetry allows doing the angular inte- 
grations, and getting rid of the 77i-quantum numbers, 
thereby making the Estates (2Z-t- l)-fold degenerate. The 
expansion yields the following Hamiltonian 

HmF = y^ (2/ + 1) {Snl - ^J.a) a\i„anlcy 
n,l.{7 

El I 

n,n' ,l,(T 

+ E^™'«in«l'U+H.c. (10) 



Here the single particle energies are £„; = 
Hu {2n + 1 + 3/2), and the Hartree interaction is 
described by the elements 






dr r'^Rniir)n„{r)Rn'i{r), 



F' , 



^„, = / drr'R^iir)A{r)R„,iir). 





We note that the a dependence of the Hartree term is 
due to the population imbalance which implies that the 
corrections are different for the two components. The 
density of a atoms is 



ncr(r) 



^ 2/ + 1 
^ 47r 

n.n'l 



Rni {r)Rn'i{r) {ali^an'ia), (11) 



and the order parameter is 



2^ + 1 



A(r) = t/(r) ^ _^i?„,(r)i?„,,(r)(al,^al,,,). (12) 



n,n' ,1 



The additional factor (2Z-I- 1) comes from the degeneracy 
of the Z-states, obtained by the summation over the m 
states. The renormalized interaction strength U{r) will 
be described later. 

Truncating the Hilbert space by keeping only the states 
with single-particle energies Eni < Ec allows writing the 
Hamiltonian in a matrix form that can be diagonalized. 
Noticing that the Hamiltonian commutes for different 
(Z)-quantum numbers makes it possible to diagonalize 
each (/)-matrix separately, simplifying the problem sig- 
nificantly. For a given (Z)-quantum numbers the Hamil- 
tonian reads 



H, 



(0 

MF 



/ -In \ 
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/ flon \ 




m' 




\ o-mi / 




V «k/i / 



(13) 



where Ni = [E/inuj) - I - 3/2] /2 is the /-specific cutoff 
(yielding the correct energy cutoff) and M' is a 2{Ni + 
1) X 2(A'';-f l)-dimensional orthogonal matrix. Notice that 
the I states have been turned into holes by switching the 
order of at and ai operators as suggested by the ordinary 
Bogoliubov transformation. The matrix M' is now 



/ £o; — Mt + ^"^ooi 



M' 



Up 



F, 



00 



^ONi 



^•^ONii 



^^N,l -MT +UJN,N,i 



^N,0 






FL 



^N.O 



-£o/ + Mi - C^^oT 



-C/J; 



OTViT 



^ON, 






(14) 



-ENil + f^l- Uf^^j^^^ J 
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(16) 
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Diagonalizing this Hamiltonian corresponds to the Bo- 
gohubov transformation, which yields the eigenener- 
gies Eji and the corresponding quasiparticle eigenstates 
(2(7V/ + l)-dimensional vectors) Wj„. The indices j, n 
both run from to 2Ni + 1. In the basis of these quasi- 
particle states, the Hamiltonian is 

27V! + 1 

H^F = E ^oi^W- (15) 

0=0 

In this basis, the density of atoms in t state is 

I j n,n'—0 

where the Fermi distribution nplE) = 1/(1 
Likewise, the density of atoms in J, state is 

"iW = E^l^E E Rni{r)Rn'i{r) 

I j n.n'—O 

X M^j,„+Ar, + iM^j,„' + Ar, + inF(-£^j7), 

(17) 

where the sign in the Fermi function is changed because 
of the I component (the last iV; + 1 components) of the 
eigenstates correspond to holes as compared to the par- 
ticles in the t components (the first Ni + 1 components) . 
The order parameter is given by 

^« = E^E E Rni{r)RMr) 

I j n,n'—0 

X l^j„W^j,„,+Ar, + i(l + 27^F(i?,,))- (18) 

The total number of atoms in the two components can 
be obtained by integrating over the densities. However, 
numerically it is faster to calculate them directly from 
the eigenstates using the equations 

2Ni+l 

N^ = E(2^ + 1) E E K^WjnMEji) (19) 



n=0 



and 



A^i 



2Ni + l 



E(2^+1)E E Kn+N,+iWl,,+N,+inF{-E,i). 

I j n=0 

(20) 



The gap and the number equations I|18ll9l20p are solved 
iteratively. We have made calculations where the Hartree 
fields are included and compared them to the case where 
Hartree fields are neglected by setting J^^'a — ^ lo'" ^11 
n, n', I . Note that close to the Feshbach resonance, the 
Hartree fields become formally infinite. In this extreme 
limit, we have limited the Hartree interaction strength U 



to 



(fepOs)" 



< /?, where P « 0.5. When comparing the 

results with and without Hartree fields, there is no differ- 
ence in the qualitative features such as the order param- 
eter oscillations and over-all shape of the gap and density 
profiles, the Hartree fields cause only minor corrections to 
gap and density profiles (effectively compressing the gas 
slightly) . Numerically, neglecting the Hartree fields gives 
a tremendous speedup in the numerical solution because 
it decouples the density and gap profiles. In the results 
presented in this article, the Hartree fields are neglected. 
We have made the calculations at zero temperature and 
present here only results at T = 0. We have checked that 
the BdG results do not change by using a finite but very 
small temperature T — 0.001 Tp. 

As the density and gap profiles are decoupled, it is 
straightforward to solve the gap equation for given chem- 
ical potentials. However, since we want to keep the num- 
ber of atoms fixed, the total procedure will require opti- 
mizing also the chemical potentials, so that the number 
equations are satisfied. The subsequent iteration proce- 
dure can be performed in several different ways. We have 
found that a very efficient procedure is to solve the chem- 
ical potentials (by solving the number equations (|19l20p ) 
for each trial gap profile A(r). The trial gap profile A(r) 
is then used for solving the new trial gap profile A'(r) us- 
ing the gap equation ifTSJ) with the new obtained chem- 
ical potential. The chemical potentials therefore keep 
changing between the iteration steps of the gap profile. 
On the other hand, the numbers of atoms stay fixed in 
the iteration process. The initial guesses for the profiles 
needed for the iteration procedure are obtained by using 
the chemical potentials and the gap profile obtained for 
a slightly lower polarization and, eventually, by the solu- 
tion obtained for unpolarized gas. Solution of the profiles 
for a given polarization P requires therefore solving the 
profiles for all lower polarizations. The validity of the fi- 
nal solution A(r), ^i, /J.^ has been checked for several val- 
ues of polarization and interaction strength by perturb- 
ing the final solution and using the perturbed profiles 
A'(r),/X|,/Lt| as initial guesses for the profiles. Usually 
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Figure 4: (color online) Gap and density profiles on the BCS 
side of the resonance (fcpas)"^ = —0.5 for P = 0.50. 



Figure 5: (color online) Typical profiles at resonance 
(kFas)'^ = 0.0 for P = 0.9. 



the disturbed profiles have converged either into the final 
profiles A(r) or into the normal state A(r) = 0. Only in 
the regime where our calculations predict the superfluid 
core polarization have our solutions sometimes strayed 
into a superfiuid solution with unpolarized core. In such 
cases, the final solution has been picked by choosing the 
solution with the lower total energy, given by 



27V, + 1 

E= ^ (2/ + 1) 



n j 

1 



U{r) 



dr r^|A(r 



(21) 



where the first (constant) term comes from switching 
to use the hole states in the I component (and this is 
countered by having half of the eigenenergies Eji neg- 
ative). The last term follows from the mean-field term 
(^i(r)^!(r))(^l(r)^l(r)) which was dropped from the 
mean-field Hamiltonian because it is a number, not op- 
erator. 

The summation in the gap equation l|18p is divergent 
without the energy cutoff Ec . This is a well known phe- 
nomenon, following from the incapability of the contact 
interaction potential to describe properly high energy 
behavior. Several different regularization schemes have 
been proposed [^,J4Q, |41|, and here we apply the one 
suggested in Ref. [4l|. This implies the following form 
for the renormalized coupling 



1 



U{r) 



1 
U 



1 

2^ 



fce(r)-fcO(0) 



■In 



fcc(r) 



where the momentum cutoff kc{r) = {2Nc + 3 — r^)^/^ 
and the local Fermi momentum fcp(r) = (/i| + fJ-i — r^)^'^ 
for the imbalanced case compared to Eq. (14) and Eq. 
(18) of Ref. [m, respectively. On the BCS side of the 
resonance for (fepOs) < we have used as the cutoff 
Ec = 200 huj and on the EEC side for {kpasy^ > the 



cutoff Ec = 240 hu! was used. We have tested the re- 
maining cutoff dependence by using the cutoff 300 hcu. 
Depending on the number of atoms, the cutoff depen- 
dence in the gap profiles was at most 2%. 



III. RESULTS 

The results show two features that we will discuss in 
detail below: 1) For small and intermediate polarizations, 
FFLO-type oscillations in the superfiuid-normal interface 
at the edge of the trap, 2) for large polarization, a polar- 
ized superfiuid that extends through the whole trapped 
gas. We study the behavior of these features, especially 
1), throughout the BCS-BEC crossover and when the sys- 
tem size, i.e. the atom number, is varied. We also calcu- 
late the condensate fraction to make a connection to ex- 
periments. Finally, we analyze which harmonic oscillator 
states are involved in pairing for different polarizations 
and discuss the connection and differences to FFLO-state 
in a homogeneous space. 



A. FFLO-type oscillations 

Typical density and gap profiles at T = 0, for small 
polarization, are shown in Figs. [Hand [21 both for LDA 
and BdG. Comparison of the profiles gives the follow- 
ing general picture: BdG predicts a) SF (equal densities) 
core, b) PS (FFLO-like oscillations) shell, c) normal state 
shell of the majority component Ni. LDA predicts a) SF 
(equal densities) core, b) normal state shell; the absence 
of the PS shell in LDA is refiected in the discontinuity of 
the density and gap profiles at the SF-N phase boundary. 
At finite temperatures, also LDA shows a polarized shell 
and the boundary becomes continuous, only showing a 
kink (Fig. [H. As shown by Figs. [iJH for small polariza- 
tion the PS given by BdG calculations is a narrow shell 
and can be understood as a boundary effect. However, 
in the following we show that, for large polarization, the 
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Figure 6: (color online) Gap and density profiles on the BEG 
side of the resonance (fcrois)"^ = 0.5 for P=0.95. 



FFLO-features extend to the center of the trap as well. 
A seen from Fig. [T]and Fig. [3] the general results from 
LDA and BdG calculation agree fairly well and the over- 
all agreement for the gap and density profiles becomes 
better for increasing particle number. The incapability 
of LDA to correctly describe the short scale behavior (ex- 
plicitly excluded by LDA) leads to the unphysical ap- 
pearance of a discontinous order parameter solution and 
can be viewed as a breakdown of LDA near the interface 
much in the same way as LDA generally breaks near the 
edge of a balanced condensate. The breakdown of LDA is 
however more pronounced in the imbalanced case due to 
presence of the normal shell surrounding the condensate 
which leads to separation of the size of the condensate 
and the size of the cloud. 



1. Dependence on the interaction strength - BCS-BEC 
crossover 

We present here results for three different cases: (1) 
BCS side of the crossover, (fcpas)^^ = —0.50, (2) unitar- 
ity, that is, (fepa^)"^ = 0, (3) BEG side (fcpas)"! = 0.50. 
We focus on the behavior at strong interactions since 
the cases (l)-(3) can all be considered to be at the uni- 
tarity regime (if it is defined \{k-pas)~^\ < 1 ) although 
representing different sides of the crossover. We do not 
consider the extreme BCS limit of weak interactions be- 
cause the features and trends observed in the unitarity 
limit are also expected to appear at weaker coupling. 

For attractive interactions, (fcpa)"^ < 0, the typical 
density and order parameter profiles looks as shown in 
Fig. |4]for polarization P = 0.50. In agreement with ear- 
lier studies using the same approach [ill, [l2, [2l|, [SJ] , the 
solution reveals an unpolarized BCS-type region at the 
center of the trap and a polarized shell with oscillating 
order parameter. The oscillations rapidly dampen when 
the density of the minority component drops. 

We have studied the presence of the oscillations around 
the unitarity region, and noticed that the critical polar- 



Figure 7: (color online) Gap profiles at {kpag) ^ — —0.5 for 
several values iV = 9000, 18000, 100000 of the total number 
of particles at a high and fixed polarization P = 0.74. For 
increasing particle number the transition from an almost con- 
stant value of the gap in the BCS core to the oscillating gap 
at the edge becomes sharper. For increasing particle number 
the gap FFLO-like oscillations at the edge the cloud display 
shorter wavelength and a slight increase in amplitude. 



ization for the appearance of the order parameter oscilla- 
tions increases with stronger interactions. At unitarity 
((fcpa)^^ = 0.0), the order parameter oscillations did 
not appear until polarization P = 0.70, compared to 
P = 0.10 at (fcpa)"^ = —0.50. The disappearance of 
the oscillations follows from the enhanced stability of the 
BCS-type pairing at stronger interactions, making the 
increased distortion of the minority component density 
profile favorable over the reduced order parameter due 
to the polarization. The result at unitarity is shown in 
Fig. [3 

On the BEC side, the oscillations do not appear. We 
have made calculations for several parameters on the 
BEC side of the resonance and have not found any FFLO- 
type oscillations in the order parameter on the BEC side. 
This is consistent with the observation discussed above 
that the critical polarization for the emergence of the 
nodes in the order parameter increases with increasing 
interaction strength. The same qualitative results on the 
BCS-BEC crossover were recently obtained in [42| using 
a hybrid BdG-LDA scheme in a cylindrically symmetric 
geometry. Typical density and gap profiles on the BEC 
side are shown in Fig. [6l 



2. Dependence on the system size - atom number N 

We have studied the dependence of the results on the 
number of atoms by performing the BdG calculations for 
the total atom numbers 9000, 10000, 18000, 20000 and 
100000. Note that this means substantial increase in the 
system size compared to our earlier work [2l| where the 
atom number 10000 was used. 

The order parameter oscillations for different numbers 
of atoms are shown in Fig. \7\ In Fermi units (i.e. in the 
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Figure 8: (color online) Exponential fits for the envelopes of 
the order parameter oscillations for 9000 and 100000 atoms. 



cloud scale Rtf), the wavelength of the oscillations be- 
comes shorter with the increased atom number but, on 
the other hand, more nodes appear. The scaling of all 
these factors is complicated, but insight can be obtained 
by fitting an exponential function ~ e~^^^ to the enve- 
lope of the order parameter oscillations as shown in Fig. 
[HI The exponential decay gives a very good description 
for the damping of the oscillations with increasing dis- 
tance from the trap center r. The penetration depths £, 
obtained from the fitting indeed slowly decrease in Fermi 
units when the atom number TV is increased, the scaling 
being roughly N^^^. Since Rtf oc 7V^/^, one may an- 
ticipate that the decrease of the penetration depth only 
occurs in relative scale, not absolute. Indeed, in a micro- 
scopic length scale (harmonic oscillator units aosc) the 
penetration depth remains constant, being 0.54 Aosc for 
100000 atoms and 0.56 Cosc for 9000 atoms. These obser- 
vations confirm that the order parameter oscillations are 
not a finite size effect but rather an interface effect for the 
superfluid-normal interface. From the shape of the order 
parameter profile it can be seen that that the interface 
formed becomes sharper for increasing particle number 
and we therefore suggest that the interface forming at 
the edge is an analog of the proximity effects appearing 
in superfiuid-normal and superfiuid-ferromagnetic junc- 
tions. In the latter case it has been shown that the inter- 
play between BCS superconductivity and ferromagnetic 
order (in our case the polarization Sn) gives rise to an os- 
cillating order parameter which decays anomalously slow 
(over several oscillations) on the ferromagnet side and 
with a characteristic length scale that is independent of 
the properties of the superconductor [43, H, [43|. An 
interesting question is the dependence of the penetra- 
tion depth on the interaction strength and polarization. 
However, we have not been able to pursue this question 
in more detail here and leave it for future work. 

In this context, we believe a few remarks on the work 
[4a | are in order. It is argued in [46|, based on a hybrid 
BdG-LDA scheme (no analysis of the penetration depths 
as presented above was done in i4Q]), that the oscillations 



Figure 9: Density and order parameter profiles at P = 0.762 
and (fepfls)"^ = —0.50 for A*' = 9000 show the superfluid with 
core polarization, together with a small dip in the density 
difference at the center of the trap. 



of the order parameter vanish for sufficiently large par- 
ticle number and are thus interpreted as a finite size ef- 
fect. It is argued that as the cutoff energy Ec introduced 
to separate the BdG and the LDA scales in the hybrid 
scheme goes to zero, LDA should be recovered. This is of 
course self-evident due to the construction of the hybrid 
algorithm but reducing the cutoff energy without con- 
cern may lead to missing some features of the system. 
The fact that a sharp interface with increasingly micro- 
scopic features (oscillations with short period on the scale 
of the cloud size) builds up for increasing number of par- 
ticles requires Ec to be increased significantly in order to 
resolve the short length scale features. Consistent with 
our results. Fig. 4 in [49| shows that the order parame- 
ter oscillates with a shorter period, and at the same time 
the amplitude of the order parameter slowly increases, 
for increasing particle number. The recent results for a 
hybrid BdG-LDA scheme [42| and the scaling analysis 
therein agree well with the results given by our full BdG 
calculations. 



B. Polarized superfluid for large imbalances 

When the polarization P exceeds some (large) critical 
value, the FFLO-type oscillations discussed above reach 
the center of the trap. At this point the gas becomes 
polarized also at the center of the trap and the order 
parameter drops to roughly half of its unpolarized BCS 
value. This realizes a superfiuid with finite polarization 
throughout the system. Fig. [9] shows the density and 
order parameter profiles for P = 0.762. For such a core 
polarized system, the concept of an interface effect is in- 
triguing as there is no BCS-type superfluid core present. 

The density difference for such a polarized superfluid 
shows an interesting feature: a small dip in the center of 
the trap (i.e. the density difference is smaller than in the 
surrounding area, but still non-zero). At zero tempera- 
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ture, this feature only appears in our BdG calculations 
whereas LDA, which fails to predict a polarized super- 
fluid at T = 0, does not lead to such a dip in the density 
difference. In contrast, LDA calculations at finite tem- 
perature produce such a feature in connection with the 
finite temperature BP phase. Recent Monte Carlo stud- 
ies of the trapped Fermi gas have shown that, for the 
strongly interacting normal state, the density difference 
increases monotonously towards the center of the trap 
[43|. Therefore one may argue that the dip is associated 
with a polarized superfiuid: at T = it is FFLO-type, at 
temperatures that are finite but clearly below the critical 
temperature it is either FFLO or BP. At higher tempera- 
tures, pseudogap effects may contribute to such features 
as well |2g|. Such a dip can be seen in the experimental 
results of |l7|. 

The interaction strength dependence of such a polar- 
ized superfiuid with oscillating order parameter is similar 
to what was already discussed above. Since there are no 
order parameter oscillations on the BEC side, we have 
also not seen any core polarized superfiuid either. Of 
course, there does exist a different kind of core polariza- 
tion on the BEC side: coexistence of a molecular conden- 
sate and free excess fermions, but that is not associated 
with oscillating order parameter. 

The parameter window for the polarized superfiuid 
shrinks with increasing atom number N. For 9000 atoms 
the window is 0.746 < P < 0.784, whereas for 18000 
atoms it is 0.746 < P < 0.774. Since the convergence 
of the calculations near this critical window is slow, we 
have not been able to determine the corresponding win- 
dow for 100000 and have not systematically analyzed how 
the window scales with particle number. For the order 
parameter oscillations, we have shown in section III. A. 2. 
that their absolute length scale stays unchanged for large 
particle numbers. Based on the present data, we cannot 
conclude with similar confidence whether the parameter 
window for the existence of the core polarized superfiuid 
becomes negligible or not for large condensates. How- 
ever, we would like to emphasize that the core polarized 
superfiuid is not due to a trivial finite size effect, i.e. not 
originating from having discrete oscillator states in the 
system description. As seen in Refs. [3J, |4l||, such fi- 
nite size effects manifest as a narrow dip in the density 
and order parameter profiles at the center of the trap. 
However, these effects vanish when the number of atoms 
increases or when the interaction strength is increased. 
Because of the stronger interactions, we do not see these 
finite size effects even at 9000 atoms. Note that such a 
dip originating from finite size effects is completely differ- 
ent from the dip in the density difference discussed above 
as a signature of the core polarized superfiuid. 



C. Condensate fraction 

To make connection to experiments [l3| where conden- 
sate fractions are measured, we calculate here the con- 
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Figure 10: The condensate fraction, density difference, and 
the order parameter as a function of polarization P obtained 
for 9000 atoms. The interaction strength is (kpa)^^ — —0.50. 



densate fraction for an imbalanced gas. In the case of a 
balanced Fermi gas the condensate fraction is defined in 
|4a . l49l | and can be viewed as a measure of the number 
of condensed pairs which in the extreme BEC limit and 
at zero temperature is just N/2 = No-. For the imbal- 
anced gas, the corresponding number of molecules in the 
asymptotic BEC limit is the number of minority atoms 
N^. It is therefore natural to consider the following nor- 
malized condensate fraction 



No/N^ 



±Jd'r^d'r2\{^^{r,)^i{r,))\' 



(22) 



for the condensate fraction in the case of N^ < Np Figs. 
[To] and [TT] show the density difference, the order param- 
eter at the center of the trap, and the imbalanced con- 
densate fraction from Eq. l(22|) as function of the polar- 
ization for (fepfls)"^ = —0.50. The critical polarization 
is roughly Pc = 0.78 and the transition is refiected in a 
rapid increase in the density difference in the center of 
the trap. Also the order parameter in the center of the 
trap drops rapidly. However, as the figures show, there 
exists a narrow but finite polarization region (in this case 
for 0.75 < P < 0.78) where both the density difference 
and the gap have non-negligible values in the center of 
the trap. All data points are for fully converged order pa- 
rameter and density profiles, satisfying the gap equation 
with precision of 10~^. In addition, several points in the 
plot have been calculated with different initial conditions 
for the iteration. 

The figures are in good qualitative agreement with the 
experimental condensate fractions [17| , showing the sud- 
den onset of the density difference at the center of the 
trap when the condensate fraction drops to zero. The 
core polarized superfiuid manifests itself as a weak re- 
vival of the condensate fraction as shown in Fig. [121 In 
other words, a finite condensate fraction co-exists with a 
finite central density difference. Although the qualitative 
agreement with the experimental results in [l3| is good, 
higher experimental accuracy as well as extending our 
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Figure 11: The condensate fraction, density difference, and 
the order parameter as a function of polarization P obtained 
for 18000 atoms. The interaction strength is (fcpa)" — 
-0.50. 



Figure 12: The core polarized superfluid region for 9000 atoms 
(fcpas)"^ = -0.50. 



IV. CONCLUSIONS AND DISCUSSION 



calculations to finite temperatures is needed for quanti- 
tative comparison, especially regarding the small window 
for the core polarized superfluid. Since the condensate 
fraction is small in the interesting transition region, bet- 
ter signatures of the polarized FFLO-type superfluid may 
be provided by the central gap, and by the dip in the cen- 
tral density difference as discussed above, provided that 
temperature dependence is also carefully investigated. 



D. Contribution of different harmonic oscillator 
states in pairing 



We studied the origin of the oscillations in the BdG 
order parameter by considering which trap shells (quan- 
tum number n) are involved in pairing in the center of 
the trap, i.e. for quantum number ^ = 0. As shown in 
Fig. [131 for zero polarization P = 0, the pairing involves 
0-2 neighboring shells, i.e. is peaked at n — n' = 0, with 
considerable weight until n—n' = 2 due to strong interac- 
tions. For small polarization P — 0.1 the peak is shifted, 
but the weight is still very much on the same n — n' as for 
P = 0, which results to only minor modiflcation of the 
order parameter proflle. However, for large polarization 
P ~ 0.75, when the oscillations of the order parameter 
appear also in the center of the trap, the pairing is peaked 
at nonzero n — n'^ 11. Therefore, for large polarizations 
the state clearly resembles the FFLO state where pair- 
ing is predominantly between momentum states k — k' 
determined by the mismatch of the Fermi surfaces (in 
our example P ~ 0.75 means a mismatch of the Fermi 
surfaces of about n — n' ~ 11). We have found that this 
behavior is not due to small particle number, in contrast, 
it becomes more clear for larger particle numbers. 



In summary, we have shown that FFLO-type oscil- 
lations appear in a trapped polarized Fermi-gas and, 
for large polarization, extend to the center of the trap 
thereby realizing a non-BCS superfluid at zero temper- 
ature. These results were obtained by BdG calculations 
using harmonic oscillator eigenstates and particle num- 
bers up to 100000. We have also made a comparison to 
the results given by local density approximation. 

The FFLO-type oscillations of the order parameter ap- 
pear on the BCS side of the BCS-BEC crossover. When 
the interaction strength is increased, the polarization re- 
quired for the existence of the oscillations grows. On the 
BEG side, no oscillating order parameter was found. The 
scaling analysis for the penetration depth of the oscilla- 
tions shows that their characteristic length scale stays 
constant when the particle number is increased. There- 
fore the characteristic length scale relative to the atom 
cloud scale (given by Rtf) decreases very slowly, ex 7V^/^. 
The features presented here should thus be observable 
even for condensate sizes corresponding to present-day 
experiments. RF-spectroscopy was proposed in our ear- 
lier work [2l| for detecting the gapless excitations related 
to the nodes of the order parameter as a signature of the 
FFLO-type oscillations. The flrst experiments using RF- 
spectroscopy in investigating imbalanced gases were re- 
cently done [50|. Moreover, the dip in the central density 
difference could provide a signature of the core polarized 
superfluid, as well as the simultaneous measurement of 
the gap and the density difference in the center of the 
trap. 

The results presented here have an interesting connec- 
tion to superconductor-ferromagnet interface effects. For 
future work, it is fascinating to think about the free- 
dom that the ultracold gases offer in terms of designable 
trapping geometries and other parameters: one should 
be able to systematically study this kind of effects from 
the limit of having large superfluid and polarized normal 
state ("ferromagnet") regions and a small interface, to the 
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limit where the interface, superfluid, and normal state are 
all of comparable size and novel mesoscopic phenomena 
may be found. 
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Figure 13: (color online) The gap at the center of the trap (the 
part of the gap for the quantum number / = ) as function 
of the difference in radial quantum number An for different 
values of the polarization and for A'' — 18000. For increasing 
polarization the mainly intra-shell and nearest-neighbor-shell 
pairing indicated by the central peak (at P = ) diminishes 
and a secondary peak appears for large An which indicate 
the increasing importance of inter-shell pairing between shell 
states having an energy separation of order determined by the 
mismatch of Fermi energies. 
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